library(MSstatsTMT)
library(tidyverse)

# batch1 -------------
# take 1-2min to import ~400M csv
batch1 <- OpenMStoMSstatsTMTFormat('~/append-ssd/nextflowing/quantms_xy_b1all_e2/msstatsconverter/xiangya.b1all.sdrf_openms_design_msstats_in.csv',use_log_file = F)

# take 17-21min to normalize???
ba1ms <- batch1 |>
  proteinSummarization(use_log_file = F)

ba1ms_prot <- ba1ms$ProteinLevelData

ba1ms_prot |>
  write_csv('mission/FPP/xiangya/sle_ms_batch1_e2_msstat.csv')

ba1ms_prot_mtx <- ba1ms_prot |>
  filter(!str_detect(Protein, ';')) |>
  as_tibble() |>
  separate_wider_delim(Condition, '|',
                       names = c('sp','ts','age','sex','dis','foo','bar')) |>
  unite('sample',c('dis','sex','age')) |>
  select(Protein,Abundance,sample) |>
  pivot_wider(names_from = sample, values_from = Abundance)

# annotate with gene symbol
b1_symbol <- ba1ms_prot_mtx |>
  separate_wider_delim(Protein, '|', names = c('type','acc','name')) |>
  pull(acc) |>
  clusterProfiler::bitr(fromType = 'UNIPROT', toType = 'SYMBOL', OrgDb = 'org.Hs.eg.db')

# save wide matrix to save space
ba1ms_prot_mtx |>
  separate_wider_delim(Protein, '|', names = c('type','UNIPROT','name')) |>
  right_join(b1_symbol) |>
  relocate(SYMBOL) |>
  select(-c(type,name)) |>
  write_csv('mission/FPP/xiangya/sle_ms_batch1_e2_wide_mtx.csv')

ba1ms_prot_mtx |>
  column_to_rownames('Protein') |>
  as.matrix() |>
  IQR(na.rm = T)

# basic qc plot
ba1ms |> dataProcessPlotsTMT(type = 'QCPlot',
                             which.Protein = 'allonly',
                             address = FALSE)

# Check the conditions in the protein level data
levels(ba1ms$ProteinLevelData$Condition)
# compare condition disease vs NC ----------
comparison <- levels(ba1ms$ProteinLevelData$Condition) |>
  as.character() |>
  str_detect('NC') |>
  if_else(-1, 1) |>
  matrix(nrow = 1)

colnames(comparison) <- unique(ba1ms$ProteinLevelData$Condition)
rownames(comparison) <- 'sce-hc'

# moderate test will throw error
test.pairwise <- groupComparisonTMT(ba1ms,
                                    contrast.matrix = comparison,
                                    use_log_file = F)

test.pairwise$ComparisonResult |>
  filter(adj.pvalue < .1)

# batch3 -------------
# take 1-2min to import ~400M csv
batch3 <- OpenMStoMSstatsTMTFormat('~/append-ssd/nextflowing/quantms_xy_b3all_e2/msstatsconverter/xiangya.b3all.sdrf_openms_design_msstats_in.csv',use_log_file = F)

# take 4min to normalize
ba3ms <- batch3 |>
  proteinSummarization(use_log_file = F)

ba3ms_prot <- ba3ms$ProteinLevelData

ba3ms_prot |>
  write_csv('mission/FPP/xiangya/sle_ms_batch3_e2_msstat.csv')

ba3ms_prot_mtx <- ba3ms_prot |>
  filter(!str_detect(Protein, ';')) |>
  as_tibble() |>
  separate_wider_delim(Condition, '|',
                       names = c('sp','ts','age','sex','dis','sample','bar')) |>
  select(Protein,Abundance,sample) |>
  pivot_wider(names_from = sample, values_from = Abundance)

# annotate with gene symbol
b3_symbol <- ba3ms_prot_mtx |>
  separate_wider_delim(Protein, '|', names = c('type','acc','name')) |>
  pull(acc) |>
  clusterProfiler::bitr(fromType = 'UNIPROT', toType = 'SYMBOL', OrgDb = 'org.Hs.eg.db')

# save wide matrix to save space
ba3ms_prot_mtx |>
  separate_wider_delim(Protein, '|', names = c('type','UNIPROT','name')) |>
  right_join(b3_symbol) |>
  relocate(SYMBOL) |>
  select(-c(type,name)) |>
  write_csv('mission/FPP/xiangya/sle_ms_batch3_e2_wide_mtx.csv')

# try merge batch1 & 3 -------------
t1 <- as_tibble(batch1)
t3 <- as_tibble(batch3)
t3$Mixture <- '3'
t3$Run <- '3_1_1'

## take 20 min to normalize
b13ms <- bind_rows(t1, t3) |>
  proteinSummarization(use_log_file = F)

b13ms |> dataProcessPlotsTMT(type = 'QCPlot',
                             which.Protein = 'allonly',
                             address = FALSE)

b13ms_prot <- b13ms$ProteinLevelData

b13ms_prot |> head()

b13ms_prot_mtx <- b13ms_prot |>
  filter(!str_detect(Protein, ';')) |>
  as_tibble() |>
  separate_wider_delim(Condition, '|',
                       names = c('sp','ts','age','sex','dis','foo','bar')) |>
  mutate(sample = case_when(Mixture == 1 ~ str_c(dis, sex, age, sep = '_'),
                            .default = foo)) |>
  select(Protein,Abundance,sample) |>
  pivot_wider(names_from = sample, values_from = Abundance)

# clearer qc boxplot
b13ms_prot |>
  filter(!str_detect(Protein, ';')) |>
  as_tibble() |>
  separate_wider_delim(Condition, '|',
                       names = c('sp','ts','age','sex','dis','foo','bar')) |>
  mutate(sample = case_when(Mixture == 1 ~ str_c(dis, sex, age, sep = '_'),
                            .default = foo)) |>
  select(Mixture,Protein,Abundance,sample) |>
  ggplot(aes(sample, Abundance)) +
  geom_boxplot() +
  coord_flip() +
  facet_wrap(~Mixture, scales = 'free')

# annotate with gene symbol
b13_symbol <- b13ms_prot_mtx |>
  separate_wider_delim(Protein, '|', names = c('type','acc','name')) |>
  pull(acc) |>
  clusterProfiler::bitr(fromType = 'UNIPROT', toType = 'SYMBOL', OrgDb = 'org.Hs.eg.db')

# save wide matrix to save space
b13ms_prot_mtx |>
  separate_wider_delim(Protein, '|', names = c('type','UNIPROT','name')) |>
  right_join(b13_symbol) |>
  relocate(SYMBOL) |>
  select(-c(type,name)) |>
  write_csv('mission/FPP/xiangya/sle_ms_batch13_e2_wide_mtx.csv')

b13ms_prot_mtx <- read_csv('mission/FPP/xiangya/sle_ms_batch13_e2_wide_mtx.csv')

mva_b13 <- b13ms_prot_mtx |>
  filter(SYMBOL %in% kegg_mva) |>
  select(-2) |>
  pivot_longer(2:last_col(), names_to = 'sample', values_to = 'abundance') |>
  mutate(group = str_extract(sample, '.+LE|NC'))

mva_b13 |>
  filter(group %in% c('SLE','NC')) |>
  ggplot(aes(group, abundance, color = group)) +
  geom_boxplot() +
  facet_wrap(~SYMBOL, scales = 'free_y') +
  scale_color_manual(values = c('blue','red')) +
  #scale_y_continuous(expand = expansion(mult = c(0,.2))) +
  theme_pubr()
  #stat_compare_means(method = 't.test', color = 'black')
